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^ ABSTRACT 

^ We numerically study the precessing disk model for superhump in the SU UMa subclass 

Q of cataclysmic variables, using a two dimensional SPH code specifically designed for thin 

disk problems. 



On 



Two disk simulations for a binary with mass ratio q = ^ (similar to OY Car) 



00 are performed, in order to investigate the Lubow (1991 a,b) tidal resonance instability 

mechanism. In the first calculation, a disk evolves under steady mass transfer from Li. 
1-H In the second simulation, mass is added in Keplerian orbit to the inner disk. The two 

CO disks follow similar evolutionary paths. However the Li stream-disk interaction is found 

to slow the disk's radial expansion and to circularise gas orbits. 

The initial eccentricity growth in our simulations is exponential at a rate slightly 
less than predicted by Lubow (1991a). We do not observe the clearing of material from 
the resonance region, via the disk's tidal response to the m = 2 component of the binary 
potential, described in Lubow (1992). Instead the m = 2 response weakens as the disk 
^ eccentricity increases. 

Oh Both disks reach an eccentric equilibrium state in which they undergo prograde 

^ precession. The rate of viscous energy dissipation in the disks has a periodic excess 

^ with a period matching the disks' rotation. The source is a large region in the outer 

^ disk. The mechanism by which the excess is produced is identified. The time taken for 

^ the periodic excess to develop is consistent with the first appearance of superhumps in 

a superoutburst. 

Key words: accretion, accretion discs - hydrodynamics - instabilities - methods: 
numerical - binaries: close - novae, cataclysmic variables 



1 INTRODUCTION 

Dwarf novae are a category of cataclysmic variables that 
exhibit semi-periodic outbursts during which the system 
brightness increases by 2 — 5 magnitudes. Cataclysmic vari- 
ables are semi-detached binaries composed of a white dwarf 
(the detached component), and a low mass secondary star 
(the contact component). For a concise introduction to cata- 
clysmic variables see Ritter (1992). In general, the secondary 
stars in cataclysmic variables are losing mass to the white 
dwarf via Roche lobe overflow. In the absence of a strong 
magnetic field, gas from the secondary forms a disk, centred 
upon the white dwarf. The source of dwarf novae outbursts 
is known to be this accretion disk. 

SU UMa systems are a subclass of dwarf novae that also 
occasionally undergo superoutbursts, which are brighter and 
longer lasting than standard outbursts. The lightcurve of an 
SU UMa in superoutburst is unique in that a component of 



the signal is modulated with a period a few percent longer 
than the binary's orbital period. This modulation is known 
as superhump. For recent superhump observations the reader 
is referred to Hessman et al. (1992), and Molnar and Kob- 
ulnicky (1992). 

The superhump signal is known to come from an ex- 
tended region in the outer part of the accretion disk (Naylor 
et al. 1987). One possible explanation, dubbed the tidal res- 
onance instability or TRI model by Lubow (1994), is that 
superhump is just the signature of an eccentric disk being 
periodically stressed by the binary's rotating tidal field. A 
slow prograde motion of the disk in the inertial frame is then 
responsible for the superhump period slightly exceeding the 
binary orbital period. In the TRI model the eccentricity is 
caused by a 3 : 1 resonance between orbits in the disk and 
the orbit of the binary itself. 

In this paper we investigate the TRI model by perform- 
ing SPH simulations of disk evolution in close binary sys- 
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terns. Firstly we review the model in more detail, contrasting 
the various simulations that were important in the model's 
development. Then we describe the adaptations needed to 
apply SPH to non-self gravitating thin disk problems. Of 
particular importance is the mechanism used to dissipate 
energy in the disk. We detail the artificial viscosity used in 
the code and provide an analytic estimate of the effective 
kinematic viscosity. This is then tested with axisymmetric 
ring spreading calculations. As a further test the result of 
a nonaxisymmetric simulation is compared against steady 
state thin disk theory. We then describe simulations of disk 
formation in the SU UMa system OY Car. 



2 PREAMBLE 

The first particle based simulations of disk formation in close 
binary systems were carried out by Lin and Pringle (1976). 
Reasoning that these accretion disks were geometrically 
thin, they restricted themselves to two dimensional simu- 
lations in the plane of the binary. Pressure forces in the disk 
were neglected. 'Viscous' forces were important and could 
not be ignored, but the details of the mechanism that pro- 
duces viscosity in these disks was not known. Thus Lin and 
Pringle employed an ad hoc technique for removing energy 
from the disk. At regular time intervals At a grid of mesh 
length I was placed over the particles in the simulation. For 
each cell the particle velocities were adjusted to place them 
in rigid rotation about their (cell) centre of mass. This was 
done in such a way as to conserve linear and angular mo- 
mentum at the cell level. 

The gravitational influence of the disk upon the binary 
was neglected, as was disk self-gravity. Thus in between 'vis- 
cous' interactions, the particles moved independently of one 
another in the potential of the binary, as test masses in the 
restricted three body problem. 

Hirose and Osaki (1990) pointed out that such a sticky 
particle scheme would have an effective kinematic viscosity 



where C was a dimensionless constant to be determined em- 
pirically. 

In this paper we use the same scalings as Lin and 
Pringle; length is scaled to the binary interstellar separation 
d, mass is scaled to the total system mass M = Ms + Mp, 
and time is scaled to the reciprocal of the binary's orbital 
angular frequency, fib~^. In these units a time interval of 
27r corresponds to one period of the binary. The geometry 
of the system is defined by the mass ratio 



Note that this definition of q agrees with most authors but 
is the inverse of the definition in Lin and Pringle. 

Although the Lin and Pringle calculations were per- 
haps hampered by insufficient resolution, they provided a 
blueprint for many subsequent simulations. Lin and Pringle 
built up their disks from zero initial mass, introducing par- 
ticles at a steady rate from the inner Lagrangian (ii) point. 
Lubow and Shu (1975) had previously calculated that gas 
overflowing the secondary Roche lobe would fall into the 



primary Roche lobe in a narrow stream and calculated its 
velocity and direction as a function of the binary mass ratio. 

As an inner boundary condition for the disk, Lin and 
Pringle simply used a hole in the domain. Any particle end- 
ing a time step less than a distance rwd from the primary 
was deleted from the calculation. Similarly any particles that 
returned to the secondary Roche lobe or that found them- 
selves a very large distance from the primary were removed 
and recorded. 

Lin and Pringle performed simulations for three mass 
ratios, q = 0.4,1.0,2.5. In each case they found the disk 
reached an equilibrium state where the rate of infall from 
Li balanced the accretion onto the primary. The total an- 
gular momentum of the disk also reached a stable value. It 
was noted that very few particles 'escaped' the white dwarf 
Roche lobe. This implied that tidal torques were the princi- 
pal means of removal of angular momentum from the outer 
edge of the disk. 

Whitehurst (1988b) implemented a fully Lagrangian 
particle method for modelling accretion disks. In contrast 
to Lin and Pringle he introduced dissipation into his code 
by allowing particles to undergo inelastic collisions with one 
another. Whitehurst also allowed for a repulsive 'pressure' 
force between particles. No grid was required to calculate 
these forces. 

The boundary conditions used by Whitehurst (particle 
removal at the white dwarf surface etc.) were similar to Lin 
and Pringle's. Whitehurst chose to make the mass input rate 
from the secondary vary with time. Pressure forces had been 
included so an equation of state was required. Whitehurst 
assumed that all the heat generated by viscous action in 
the disk was instantaneously radiated away. Thus he set the 
speed of sound to a constant value. 

Whitehurst (1988b) described a calculation for a binary 
with q = j^. The disk developed a two armed spiral struc- 
ture and was elongated perpendicular to the binary axis. 
The disk appeared stationary in the binary frame. The mean 
outer radius of this disk corresponded with Paczyhski's 
(1977) estimate. Paczyhski had assumed that gas stream- 
lines in the disk would be well approximated by a family 
of simply periodic orbits of noninteracting particles in the 
restricted three body problem. Close to the primary these 
orbits are approximately circular and concentric. At larger 
distances the orbits are more distorted and at larger radii 
again the orbits begin to intersect one another. Paczyhski 
proposed that the disk would be truncated at the last nonin- 
tersecting simply periodic orbit of a noninteracting particle 
(LNO). This is referred to as tidal truncation of the disk. 

A simulation of a system with a low mass secondary 
(q=0.15, Whitehurst 1988a) produced a more eccentric disk. 
In the binary frame the disk's semi-major axis rotated in the 
opposite sense to the motion of the disk gas with a period 
slightly exceeding the binary period. (In the inertial frame 
the disk exhibited a slow prograde precession). The energy 
dissipation in the disk had a cyclic component correspond- 
ing to the disk rotation period. Whitehurst proposed that 
superhump in SU UMa binaries was due to the same tidal 
instability seen in his simulation. The superhump period 
correlated with the rotation of the simulation disk. Also, the 
tidal instability only occurred in low mass ratio systems, and 
SU UMa systems characteristically have low mass ratios. 

Hirose and Osaki (1990) performed a series of two di- 
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mensional simulations for different mass ratios using tlie Lin 
and Pringle sticky particle method. In each calculation they 
allowed a disk to form from scratch via steady Taa,ss transfer 
from the secondary. They found that in general, disks in low 
mass ratio systems were tidally unstable whilst those in high 
mass ratio systems were stable. 

The range of mass ratios for which tidal instability oc- 
curred depended upon the resolution (and effective viscos- 
ity) of their simulations. One series of simulations with grid 
size I = O.Of (using the same scalings as above) produced 
stationary disks for mass ratios q > 0.25. A second set of 
calculations with I = 0.02 rf, and hence with four times the 
effective viscosity, gave forth tidally unstable disks in all but 
the q = 1 case. In an argument that built upon Whitehurst's 
hypothesis, Hirose and Osaki proposed that the tidal insta- 
bility was caused by a 3 : f resonance between orbits in the 
disk and the orbit of the binary itself. With an analysis of 
the orbits of noninteracting particles they found that the 
3 : f resonance only occurred inside the tidal radius (LNO) 
for systems with q ^ 0.25. Hirose and Osaki rationalised 
the differing stability results of the two simulation sets by 
explaining that in the I = 0.02 d calculations, the large effec- 
tive viscosity caused the disks to spread significantly beyond 
the tidal radius, allowing them to interact with the 3 : f res- 
onance in larger mass ratio systems. 

Lubow (f 99f a) added a good deal to the sophistication 
of the model. Rather than considering the orbits of noninter- 
acting particles, he completed an analysis of a fluid disk or 
ring's response to a tidal field. Knowing that eccentric Lind- 
blad resonances cause eccentricity growth, he identified the 
m = 3 eccentric inner Lindblad resonance as being the most 
important resonance for increasing eccentricity in accretion 
disks in close binary systems. With a nonlinear mode analy- 
sis Lubow went on to identify the feedback mechanism caus- 
ing the eccentricity to grow and then calculated a growth 
rate. 

In a set of cylindrical coordinates (r, 9) centred upon 
the primary, the disk eccentricity can be expressed as a si- 
nusoidal mode with argument 9. The tidal field due to the 
secondary star can be decomposed into a set of functions 
4>m{r)cos[m(9 — fibi)], where m is an integer. The m = 3 
component of the tidal field generates a sinusoidal response 
in the disk with argument 3^ — 3Qbi. At the Lindblad res- 
onance this response couples with any existing eccentricity 
to produce a travelling wave (argument 29 — 3Qbi). In turn, 
the interaction of the wave and the tidal response acts to 
reinforce the disk's eccentricity, which grows exponentially 
with a rate 

A = 2.08x27rQbg^r,es^^%^#^. (3) 

Md (e) 

Here rres is the radius of the Lindblad resonance, E(r) is the 
surface density at radius r, Md is the total disk mass, e(r) 
is the eccentricity at radius r, and (e) is the mass averaged 
disk eccentricity. 

In a second paper, Lubow (f 99fb) detailed several SPH 
calculations that confirmed both the mechanism and the cal- 
culated growth rate. His simulations however differed from 
those of Whitehurst, and from those of Hirose and Osaki. 
Rather than allowing a disk to build up via mass transfer 
from the secondary, Lubow prescribed an initial disk with 
particles laid down on the simple periodic orbits of noninter- 



acting particles. The mass transfer stream was completely 
absent from Lubow's calculations. Secondly, in order to view 
the newly discovered mechanism in isolation, Lubow initially 
included only the m = 3 component of the tidal field in 
his calculations. He found that when the full binary poten- 
tial was included, the eccentricity growth rate was reduced. 
He later determined (Lubow f992) that non-resonant tidal 
stresses due to the m = 2 component of the tidal field re- 
pelled material from the inner Lindblad resonance and hence 
reduced the instability's strength. Lubow expressed concern 
that the weakened instability would be incapable of produc- 
ing a superhump response sufficiently rapidly after the onset 
of superoutburst. 

The major factors contributing to the disk's preces- 
sion (a slow prograde motion with respect to the inertial 
frame that, in the TRl model, is responsible for the super- 
hump period slightly exceeding the binary period) were de- 
tailed in Lubow (f992). The m = component of the tidal 
field induces a prograde precession by shifting the epicyclic 
frequency from the Keplerian value, Qb. Secondly, Lubow 
found that pressure waves would produce a retrograde pre- 
cession. A third factor, due to wave stresses, was found to 
be capable of producing both a prograde or a retrograde 
procession. 

The analysis up till that point did not consider the 
possible effects of a sudden, dramatic increase in viscos- 
ity in the disk that would occur in the thermal instability 
model of a dwarf novae outburst. Neither did the model ac- 
count for influence of the gas stream from the secondary. 
Lubow (f994) did examine the role the stream might play 
in creating or increasing eccentricity in the disk. In partic- 
ular he considered the periodic variation in the momentum 
the gas stream would impart to a precessing elliptical disk 
and asked whether it might act to increase the eccentricity. 
He found in fact, that it would reduce the disk's eccentric- 
ity. The analysis did not consider the distortions upon disk 
streamlines that the secondary would impart. 

There is a second question to be asked with regard to 
the gas stream. We know that, in comparison to gas in the 
outer disk, the gas stream has relatively low angular momen- 
tum and their collision will reduce specific angular momen- 
tum in the outer disk. So does the stream disk collision act 
to push gas into the resonant region? We will try to answer 
this question in this paper. 

To reiterate, disagreement over the Tidal Resonance In- 
stability (TRl) model for superhump appears to be focused 
upon whether the gas is capable of reaching the 3 : f res- 
onance. If that issue is to be resolved then we must know 
more about the mechanism dissipating energy and trans- 
porting angular momentum through the disk. 

To emphasise the previous point we consider the sim- 
ulations of Heemskerk (f994). He assumed the gas to be 
inviscid, and conducted simulations on that basis using an 
Eulerian code. In calculations he performed that included 
only the m = 3 component of the secondary's potential, 
he found the eccentricity grew in agreement with Lubow's 
model. However, when Heemskerk included the full binary 
potential in his simulations, he found the disk to be trun- 
cated somewhat inward of the resonance region and that 
there was no eccentricity growth. But, apart from any nu- 
merical viscosity present in his code, Heemskerk had not 
explicitly included a mechanism for transporting angular 
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momentum through the disk. It is not surprising therefore 
that his disks should have a smaller outer dimension than 
the highly viscous disk calculations of say Hirose and Osaki. 
As an aside we remark that Heemskerk was incorrect in re- 
ferring to the calculations of Whitehurst and of Hirose and 
Osaki as being SPH simulations. The sticky particle method 
used by Lin and Pringle, and Hirose and Osaki relies upon 
a grid for determining the energy dissipation. SPH on the 
other hand is characterised by its independence of any mesh 
scheme. And whilst Whitehurst's unique particle scheme was 
also fully Lagrangian, it did not explicitly include the kernel 
interpolation scheme that is central to SPH. 



3 NUMERICAL METHOD 

The numerical simulations in this paper were completed us- 
ing Smoothed Particle Hydrodynamics (SPH), a fully La- 
grangian particle method. For a comprehuensive review of 
SPH, the reader is referred to Monaghan (1992). In this sec- 
tion we discuss the specific application of the SPH technique 
to accretion disk problems. First, an appropriate set of SPH 
equations for geometrically thin disks is obtained. Then dis- 
sipation in the SPH scheme is discussed. The standard ar- 
tificial viscosity term is modified to a form more suitable to 
the modelling of viscous shear, and a theoretical estimate 
of its equivalent kinematic viscosity is obtained. The esti- 
mate is tested in the next section. The artificial viscosity is 
related to the Shakura-Sunyaev viscosity parameterisation 
commonly used in thin disk theory. Finally we consider the 
fact that viscous time scales in an accretion disk are typi- 
cally much longer than dynamical time scales, and discuss 
the problems this presents for numerical modelling. 



3.1 SPH Equations for a Geometrically Thin Disk 

We start from the basic equation for an integral interpolant. 



Ai{r)= A{r') W{r -r',h) dr. 



(4) 



W{r, h) is an interpolating kernel, and h is the smoothing 
length. The functions A that we are interested in describe 
disk properties that have in some sense been vertically inte- 
grated to account for the disk thickness, and so are functions 
defined on the disk midplane. Accordingly we will use cylin- 
drical coordinates {r{r, 9), z) with z = on the midplane. 
An element of disk surface area dr with surface density E 
has mass 



dr 



(5) 



As usual in SPH the integral interpolant is approximated by 
a summation over a set of disordered particles. 



Asir) = V mt^W{r - rt, h), 
^ — ' Eb 

b 



(6) 



with nib being the mass of particle b. Thus the surface den- 
sity at r is given by 



Es(r) = ^ mbW{r - Vb, h). 



(7) 



To (approximately) find the volume density p on the mid- 
plane we divide E by the disk's pressure scale height, which 
is given by 

^ = i- 

We now want the SPH equivalents for the momentum and 
energy equations for thin disks. For particle a, the standard 
SPH momentum equation incorporating artificial viscosity 



dva ,Pa , Pb 

b 

+ ilh!^^^^^!i:}h±)w(Ya-Vb,h). 

Pab 



(9) 



Here Cab = 0.5(ca -|- Cb) is the average sound speed for parti- 
cles a and h. pab is similarly defined. The notation is made 
more concise by writing FF(ra — rb, h) as Wab- C ^nd /? are 
the linear and nonlinear artificial viscosity parameters re- 
spectively. 



r^^2_|.^2 Vab lab ^ U 



P-ab 



(10) 







Vab • rab > 



where Vab = Va — Vb and rab = ra — rb. is a constant with 
the dimensions of length. Typically = 0.1 h. 

We take the (r, 9) component of the fluid momentum 
equation. 



At p 



(11) 



Following Goldreich and Lynden-Bell (1965) we multiply 
Eqn 11 by /9 and integrate over z to obtain 



dvr,e , 



(12) 



Reversing the order of operations on the right hand side, 
and assuming Vr,e is independent of z, we get 



dVr 



dt 



+ -V,.e / Pdz 



^^r.eiPH) 



We rewrite the pressure gradient term, 
- Vr.e (PH) = Vr.e ( — ) + — Vr.e E. 



(13) 
(14) 

(15) 



Ignoring the artificial viscosity term for a moment, our two 
dimensional SPH momentum equation for a geometrically 
thin disk is 



^ = - 2^ + ^) VaPFab. 

b ° 

Similarly, our SPH energy equation is 



dMa 1 



mb(- 



,PaHa , PbHb 



) Vab • VaWab. 



(16) 



(17) 
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3.2 SPH Artificial Viscosity 



The standard SPH artificial viscosity expression, as seen in 
Monaghan, is composed of a linear term and a quadratic 
term. The quadratic (/?) term, which is similar to a Von 
Neumann- Richtmyer viscosity, was originally included to im- 
prove the treatment of high Mach number shocks. We are 
interested in a supersonic shearing flow for which the linear 
(^) term alone is sufficient. Thus we have set /? = 0, as did 
Meglicki et al. (f993). Traditionally a is used to denote the 
SPH linear artificial viscosity parameter. We have changed 
the notation to avoid confusion with the Shakura-Sunyaev 
parameter. 

In the standard formulation, the artificial viscosity term 
is set to zero for expanding flows. Again we realise that the 
viscous term was originally intended to improve SPH's han- 
dling of compression, and dissipation in regions of rarefac- 
tion was undesirable. Here however, with a strongly shear- 
ing flow, it makes no sense to use such a viscous switch. 
The remaining linear ^ term generates both shear and bulk 
viscosity. 

By reversing the procedure used to generate the SPH 
summation equations we can obtain a continuum limit for 
this artificial viscosity term (see e.g. Pongracic (f988), or 
Murray (f994)). When we let the number of particles N 
oo and the smoothing length A ^ we obtain 



C^hnc 
2E 



(V-(cES)-F V(cEV-v)), 



where the deformation tensor 
_ dvg dvb 

OXb OXa 



(18) 



(19) 



K is a constant determined by the interpolation kernel being 
used. In two dimensions, with the cubic spline kernel given 
in Monaghan (f992), we find 



(20) 



If we assume the density and sound speed vary on length 
scales much longer than the velocity we have 



av = ^(V"v + 2V(V.v)). 



Figure 1. Radial surface density plots at times t = 0.0,0.2,0.4 
and 0.8 for the axi symmetric ring spreading calculation. Param- 
eters used in the simulation are listed in Table 1. Eqn 32 is also 
plotted (solid lines). For this simulation, particles were laid down 
in 21 concentric rings which were •^h apart. 



(21) 



Figure 2. Maximum surface density of the simulation depicted 
in Fig. 1, as a function of time (heavy dots). Eqn 32 evaluated at 
r = tq (solid line) is plotted for comparison. 



The simulations described later in this paper have constant 
c, C, and h. In terms of the Shakura-Sunyaev parameterisa- 

_ 3. 

tion, this is like setting a oc r 2 . 



If we assume the velocity divergence is zero then we have 

(22) 



Che „2 



Thus in the continuum limit we could expect the linear ar- 
tificial viscosity term to generate a shear viscosity 



-Cch. 



(23) 



The form of Eqn 23 is similar to the Shakura Sunyaev vis- 
cosity parametrisation 



acH. 



(24) 



In the SPH case however, the length scale over which dissi- 
pation occurs is the smoothing length, rather than the pres- 
sure scale height. A disk with kinematic viscosity given by 
Eqn 23 has a viscous dissipation equivalent to an a disk with 



1 + g' 



C - 

— r 

c 



(25) 



3.3 Axisymmetric Tests 

In order to verify our analysis of the standard SPH artificial 
viscosity term we completed a series of axisymmetric ring 
spreading simulations for comparison with semianalytic so- 
lutions. Lynden-Bell and Pringle (f974) considered the evo- 



Table 1. Parameter settings used in tlie axisymmetric ring 
spreading calculations 



Parameter 


Value 


Remark 


c 


10.0 


SPH artificial viscosity parameter 


h 


0.01 


SPH smoothing length 


c 


0.02 


Sound speed 


ro 


0.85 


centre of Gaussian 


ri 


0.80 


inner edge of annulus 


r2 


0.90 


outer edge of annulus 


I 


0.025 


half width of Gaussian 


Ar 




initial interparticle spacing 
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poo , 

I Mr') exp[-4i^] /.(^)d/. (30) 

Upon trying several initial density functions, we found that 
a smoothly varying Eo gave best results. We thus chose a 
Gaussian, 



Figure 3. Maximum smface density of a second ring spreading 
simulation as a function of time (heavy dots). Parameters are 
identical to the previous simulation with the exception that Ar = 
h. Again Eqn 32 evaluated at r = ro (solid line) is plotted for 
comparison. 



lution of an axisymmetric ring with an initial density distri- 
bution 



Eo(r) 



27rro 



6{r — ro), 



(26) 



where 8 is the Dirac delta function and ro is the ring's initial 
radius. For a constant kinematic viscosity j/, and a Keplerian 
potential, the solution for t > (Eqn 2.13, Pringle 1981) is 



E(r,t) 



exp 



(1 + 



/.(-), 
4 r 



(27) 



where x = — and r = ^-^^ are dimensionless radius and 
. ° . 

time variables respectively, and Ix is a modified Bessel func- 

4 

tion of the first kind. Now for large arguments 



/2wx 

So for small r, 
E(r,t) 



X > 1. 



1 



27r 2 To x' 



■ exp 



(1 



r < 1. 



(28) 



(29) 



The density distribution at some small time r is thus ap- 
proximately a Gaussian modified by factor x~^^^. 

Flebbe et al. (1994) sought to test their SPH code, 
which incorporated a tensor form of a general Navier Stokes 
viscosity term, by simulating such a ring. Naturally with 
SPH they could not set up the initial 8 density distribu- 
tion so they started the simulation at dimensionless time 
r = 1.6 X 10~* . It is not clear how Flebbe et al. constructed 
their ring to obtain the correct density profile. In general one 
has the option of either manipulating the particle number 
density in the ring or initialising the particles with different 
masses. 

Rather than simulate the evolution of a 5 function ring, 
we chose to make use of the fact that when is constant 
(or indeed when is a function of r only) the equation for 
the evolution of an axisymmetric ring is linear in the density 
and a solution for a general initial condition can be obtained 
from the 8 function response. Thus 

E(r,t)= ^ 



Eo(r) 



exp {- {^—f^Y} Tl <T <T2 

elsewhere 



(31) 



truncated at ri and r2. ro and I are constants specifying the 
radius of maximum density and the width of the Gaussian 
curve respectively. Using Eqn 28 we obtain, for small time 
t, the approximate solution 



E(r,t) 



r3/1\/irl2i/t 



exp 



exp 



(r' 



12vt 



dr'. (32) 



Once the density is known, the radial component of the ve- 
locity is given by 



Vr(T,i) = -3f^ |ln(r2E(r,t)) 



For Eo given by Eqn 31 we have 



Vr{r, 0) 



31/ 61/ , 

1 (r ■ 

2r |2 I 



ro 



(33) 



(34) 



In order to initialise our ring we require the kinematic vis- 
cosity. In fact we simply substituted Eqn 23 for v in Eqn 34. 

In the simulations described here, the kinematic viscos- 
ity was assumed to be constant. The code used a smooth- 
ing length h that was constant in time and space. Pressure 
forces were switched off so we could study the artificial vis- 
cosity term in isolation. However, as oc c, the sound speed 
was not set to zero. Best simulation results were obtained 
when the initial particle distribution reflected the symmetry 
of the system. Particles were laid out in a series of concen- 
tric rings Ar apart, so as to give a constant number density 
throughout the ring. Particle masses were initialised to give 
the Gaussian density profile. To provide comparison with 
simulation results, Eqn 32 was numerically integrated using 
MATHEMATICA. 

In Fig. 1, simulation density profiles (heavy dots) at 
the times shown are compared with Eqn 32. The various 
parameters used in this simulation are summarised in Ta- 
ble 1. Units were chosen so that the mass of the central 
body, and the orbital angular frequency of a test particle 
at unit radius, were both equal to one. Flebbe et al. (1994) 
stated that in their simulation the radial velocity was less 
than 5% of the azimuthal velocity. In our case the figure is 
approximately 10%. Eqn 25 indicates that this simulation 
generates a dissipation comparable with a Shakura-Sunyaev 
viscosity with parameter a = 0.8. This is not to say that the 
SPH artificial viscosity is producing a dissipation equivalent 
to a Shakura-Sunyaev term. We merely include the figure 
for order of magnitude comparison. 

Densities in the simulation were calculated using Eqn 7, 
as opposed to integrating the continuity equation. As a re- 
sult the initial peak density in the simulation is approxi- 
mately 5% less than that of the Gaussian. This is because 
h is comparable to the Gaussian half width 1. At the peak 
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the density gradient changes significantly over a smoothing 
length. Improved resolution could be gained by reducing h 
or increasing I. A similar simulation with smaller h would 
require an increased number of particles, and hence more 
computer time. On the other hand the ring spreads with a 
time scale that is proportional to P so increasing the Gaus- 
sian half width would also lead to more time consuming sim- 
ulations. Fig. 1 may be compared to Fig. 2 of Flebbe et al. 
Their results show more scatter about the analytical curves, 
possibly as a result of the way in which their particles were 
set up. 

In Fig. 2 we plot the ring's peak density as a function 
of time (heavy dots). We also show E(ro, i), i.e. Eqn 32 eval- 
uated at To- For small time t this gives a good estimate of 
the ring's maximum density. At first glance the agreement 
between theory and simulation does not appear as good as 
Fig. 1 suggested. We have already discussed why the simu- 
lation peak density should be slightly reduced at small time. 
At times t > 1 the simulation density is larger than the the- 
oretical value and appears to be levelling off more rapidly. 
As the annulus spreads, the gaps between the constituent 
rings of particles get larger until finally they are comparable 
with the smoothing length and the rings lose contact with 
one another. At that point, a particle's only near neighbours 
lie on the same ring. This is a resolution problem which can 
be overcome by increasing the number of particles or alter- 
natively by incorporating a variable smoothing length. To 
confirm our explanation we repeated the simulation with an 
increased number of particles (Ar was reduced from i/t to 
Y^h). Fig. 3 shows shows that the resulting peak density 
profile agrees much better with the semianalytic curve. 

Artymowicz and Lubow (1994) used similar tests to es- 
timate the shear viscosity of their SPH code. In contrast to 
our artificial viscosity term which is added to the pressure 
term in the momentum equation, Artymowicz and Lubow 
choose a multiplicative term. The nature of the term is es- 
sentially the same as ours, differing only by a factor to 
maintain dimensional consistency. Artymowicz and Lubow 
also found it appropriate to remove the viscous switch. Pon- 
gracic took a term similar to Artymowicz and Lubow's to 
the continuum limit, in the manner described in the previ- 
ous section, and obtained an equation similar but not identi- 
cal to Eqn 18. Once the appropriate assumptions have been 
made however, such a term would also produce the shear vis- 
cosity given by Eqn 23. Without the aid of such an analysis, 
Artymowicz and Lubow obtained via dimensional arguments 
the viscosity relation 



f = gCc(/i), 



(35) 



with q being a constant of proportionality obtained empiri- 
cally from ring spreading calculations. As their code used a 
variable smoothing length, Artymowicz and Lubow used a 
locally averaged smoothing length (h) in Eqn 35. For simu- 
lations using a three dimensional code with the cubic spline 
kernel they found q = 0.1 to within 30%. This would appear 
to be more in agreement with the analysis for the two di- 
mensional rather than the three dimensional code. However 
we note that because the Artymowicz and Lubow artificial 
viscosity term multiplies the pressure term in the momen- 
tum equation, they could not switch the pressure forces off 
to consider the effects of the viscous dissipation term in iso- 
lation. The pressure causes an additional spreading of the 



Figure 4. Azimuthally averaged surface density plotted as a func- 

1000.00 (heavy 



tion of r for the 



1 simulation at t 



dots). Simulation parameters are listed in Table 2. Eqn 36 with 
r, = 0.017 = 2.35 X lO"* d^fi"^ gives best fit to the data and 
is plotted as the solid line. 



Figure 5. Azimuthally averaged rate of energy dissipation in the 
g = 1 simulation at t = 1000.00 Qj^"^ plotted on logarithmic axes 
(heavy dots). Eqn 38 with = O.OlTrf is also shown (solid line). 



ring which leads to an error in estimating q. Artymowicz 
and Lubow do not go into the details of their tests, and in 
particular it is not clear how they chose their initial radial 
velocity condition. 

We found that if, instead of using the kinematic viscos- 
ity given by Eqn 23, we used a different to determine the 
ring's initial radial velocity, the agreement between simula- 
tion and theory deteriorated. Thus we were able to confirm 
that Eqn 23 gave an accurate estimation of the viscosity, to 
within an uncertanity of approximately 10%. 



3.4 Nonaxisymmetric Tests 

As discussed in Section 2, simulations by Lin and Pringle 
(1976), Whitehurst (1988b), and Hirose and Osaki (1990) 
for high mass ratio systems had resulted in steady state 
disks that were approximately axisymmetric and stationary 
in the binary frame. In fact Hirose and Osaki used such 
a simulation to determine the viscous dissipation of their 
sticky particle code. 

The surface density in a steady state axisymmetric disk 
is given by 
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Table 2. Parameters for g = 1 disk simulation. 



Parameter 


Value 


Remark 


c 


10.0 


SPH artificial viscosity parameter 


h 


0.01 d 


SPH smoothing length 


c 


0.02 dur^ 


Sound speed 


■Winj 


0.1 


velocity of particles injected at Li 


''inj 


0.3491 rad 


angle to binary axis 






of injected particles 


At 


0.01 


interval between particle injection 




0.02 d 


white dwarf radius 



,E = ^fl-,/^, (36) 
37r V V r ; ' ^ ' 

where M is the steady state mass flux through the disk and 
r, is the radius of the disk's inner edge. The inner boundary 
condition implicit in Eqn 36 is 



Hirose and Osaki took the steady final state of their simula- 
tion and fit the density distribution to Eqn 36, using r, and 
J/ as free parameters. 

The rate at which energy is dissipated per unit area in 
a steady state axisymmetric disk 

= l^(l-.f^. (38) 

Hirose and Osaki compared energy dissipation in their sim- 
ulation with Eqn 38. In this case the single free parameter is 
r,. The dissipation comparison thus provided a check upon 
the density comparison. 

For us, the nonaxisymmetric test provides a useful in- 
termediate step between the ring spreading calculations and 
the simulations we are interested in finally performing. Full 
disk simulations typically take several tens of thousands of 
time steps and there is no direct way of measuring their ac- 
curacy or the validity of the result. With the particle number 
constantly changing and viscous forces at work, there are no 
constants of the motion that can be monitored as the simu- 
lation proceeds. In contrast the ring spreading calculations 
could be completed in a hundred or so time steps and we 
were able to follow monitor conservation of angular momen- 
tum. We briefly describe therefore a simulation in a system 
with mass ratio q = 1. 

We use similar boundary conditions to Lin and Pringle 
(1976). In particular we choose the simplest inner boundary 
condition, that is a circular hole in the domain of radius rwd 
centred upon the primary. Any particle ending a time step 
within this circle is simply dropped from the calculation. We 
record and remove any particles that return to the secondary 
or that are flung to a large radius. 

Initially, simulations were run with the pressure forces 
set to zero, the only interparticle forces being due to the ar- 
tificial viscosity term. It was found however that a number of 
tightly bound particle pairs formed. The artificial viscosity 
force between two particles always opposes the relative mo- 
tion. Thus it was possible for two particles to approach each 
other very closely and become 'stuck' to one another. The 
fraction of particles that formed such pairs was only ever 
very small (less than half a percent). But they produced 
density extremes which caused the Courant condition time 



step to become very small, and hence the simulation to be- 
come very slow. However pressure forces, which are always 
repulsive, act to prevent the formation of such pairs. Thus 
the pressure terms were included in all further simulations. 

The simulation began with no particles in the disk. 
Then, following Hirose and Osaki, single particles were in- 
jected at regular intervals At from the inner Lagrangian 
point. Simulation parameters are summarised in Table 2. 
For these values Eqn 23 predicts u = 2.5 x 10~* Qb. 

We followed the calculation to t = 1000.00 Qb"^ (159 
binary periods). The disk reached a steady state after about 
100 binary periods. 

The radial density profile for the simulation disk at 
t = 1000.00 Qb"^ is plotted in Fig 4, along with the best 
fit theoretical curve (r, = 0.017 = 2.35 x 10"* rf^ Qb). 
The q = 1 steady state profile agrees with axisymmetric 
theory out to a radius r ~ 0.24 rf. At r ~ 0.26 there is a 
peak in the simulation curve due to a gas buildup or dense 
'ridge' that runs partway around the outer edge of the disk, 
coinciding with Paczynski's (1977) last nonintersecting or- 
bit. The small peak at r ~ 0.35 is due to the compression 
of the disk by the injection stream at the disk hot spot. 
The r > 0.4 component of the curve is due mainly to the 
injection stream. 

The viscosity in the simulation agrees to within 6% with 
Eqn 23. We found that only a small range of values for (i.e. 
±5 X 10~^ Qb) gave a reasonableiit to the steady state disk 
profile. We thus have a second independent confirmation of 
Eqn 23 to within about 10%. Note that pressure forces were 
included in the simulation, which may explain some of the 
small discrepancy. 

The value of r, giving best fit is slightly less than rwd. 
This is not worrying considering our open inner boundary. 

The energy dissipation in the simulation is estimated us- 
ing the SPH energy equation. In fact, the equation accounts 
for both pressure and viscous forces so in regions of expan- 
sionary flow with little shear, ^ may be negative. Thus 
comparison of simulation results with Eqn 38 also yields in- 
formation regarding the relative importance of pressure and 
viscous forces in the simulation. The radial distribution of 
energy dissipation in the disk at t = 1000.00 Qb~^ is shown 
in Fig. 5 along with Eqn 38 (r, = 0.017 rf). There is very 
good agreement for the inner disk, revealing the predomi- 
nance of viscous forces there. 



4 RESULTS 
4.1 Simulations 

Two simulations were completed. System parameters were 
chosen to correspond with the eclipsing SU UMa system 
OY Car. Based upon the observations of Wood et al. (1989), 
a mass ratio q = -p^ was used. The binary orbital period 
T = 0.063121 days. A steady mass transfer rate from the 
secondary of 10~^MQyr~"' was assumed. According to the 
eclipse mapping observations of Rutten et al. (1992) this is 
a typical M for OY Car in outburst. An isothermal equation 
of state 

P 2 

— = c = constant (39) 
P 
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Figure 6. Particle number as a function of time for simulations 
1 and 2. Using the system parameters for OY Car given in the 
text, an SPH particle has a mass 0.4 X 10""^^ times the binary 
mass. Hence this figure can be read as a disk mass versus time 
graph. 



was used. Other calculations with different equations of state 
have been performed but are not detailed here. They include 
simulations of constant Mach number disks, and simulations 
of disks whose surface elements radiate as blackbodies. 

With the exceptions noted below, these two simulations 
were completed in a very similar fashion to the nonaxisym- 
metric test in Section 3.4. Simulation parameters are given 
by Table 2. Once more th sound speed c = 0.02 rfQb and 
viscosity p = 2.5 x 10~* rf^ Qb. Again the calculations were 
followed till t = 1000.00 Qj^^. Boundary conditions remain 
unchanged. 

The two q = simulations differed only in the details 
of the particle injection: 

Simulation 1 Particles were added singly at Li every 
At = 0.01 Qb~^ with a speed v = O.lrfQb at an angle .367 
rad to the binary axis in the direction of the binary rotation. 

Simulation 2 Particles were added to the calculation by 
placing them in a circular Keplerian orbit at the circulari- 
sation radius, Tc- This is just the radius of the narrow ring 
about the white dwarf that forms early on in Simulation 1. 



For q 



0.1781 d. Thus mass and angular momen- 



tum was added to the disk at the same rate as in the first 
simulation. 



4.2 Evolution 

In Fig. 6, the disk mass is plotted as a function of time 
for the two simulations. In both cases, the mass increases 
to a maximum before decaying slightly to a steady value. 
A sticky particle simulation of Hirose and Osaki (1990) for 
a system with q = 0.15 and = 2.02 x 10~* Qb had 
a similar mass-time diagram (their Fig. 5). For all three 
calculations the evolutionary picture is as follows. Initially 
the disk remains approximately axisymmetric as it increases 
in mass. Then the 3 : 1 resonance is encountered and the 
disk develops an eccentricity. In the binary frame the disk's 
major axis begins to rotate in a retrograde fashion with a 



period a few percent longer than the binary's (in the inertial 
frame the motion is prograde). 

In these three simulations the disk became sufficiently 
eccentric that when the disk centre of mass passed between 
the two stars, the outer edge just brushed the secondary 
star's Roche lobe. Maximum disk mass coincided with the 
onset of this cyclic mass return to the secondary. The disk 
then finds an equilibrium between mass input from Li, ac- 
cretion onto the white dwarf, and mass return to the sec- 
ondary. Note that in our simulations the ratio of mass ac- 
creted onto the white dwarf to mass returned to the sec- 
ondary was greater than 100 : 1. 

The simulation 1 disk encountered the 3 : 1 resonance at 
a much later time than did the simulation 2 disk. Clearly the 
addition of material with low specific angular momentum to 
the outer disk inhibited the disk's expansion to larger radii, 
resulting in a much larger steady state mass for simulation 
1. 

Figs 7 and 8 show the surface density at various epochs 
for simulations 1 and 2 respectively. The mass transfer 
stream's role as a 'growth inhibitor' is readily apparent. 
We draw particular attention to the relative states of the 
two disks at time t = 50.00 Qj^^ (top left frame). Disk 1 is 
confined to an almost circular region well within the tidal 
radius. There is a sharp drop in the density at the disk's 
well defined outer boundary. In contrast, disk 2 extends out 
to the tidal radius with a very shallow density gradient in 
the outer regions. 

Early on in the development of the tidal instability, two 
armed spirals are evident in both simulations (frame 4 Fig. 7 
and frame 2 Fig. 8). We identify these as the {20 - 3Qbt) 
mode that is a component of the tidal instability mechanism. 
Later in the simulations, with several modes excited, it is 
harder to identify them individually. 

4.3 Eccentricity growth rate 

To determine the strength of the W — mQb^ mode in the 
disk, hereafter referred to as the (I, m) mode, we follow 
Lubow (1991 b) and Fourier decompose the simulation par- 
ticle distributions as functions of angular position and time. 
We choose a non-rotating coordinate system (r, 9) centred 
upon the white dwarf. By definition, the total strength of 
the (I, m) mode 

Sl,m{i) — {Scos,cos,l,m{i) ~t" '5'cos ,sin ,^ ,m (^ ) 
2 , o 



+ 

in ,sin ,^ ,m (t ) ~t~ '5'sin ,cos ,^ ,m (^ 

The component modes are defined as follows 

2 



(40) 



t + 2jr 



E 



N 

p=l 



7rN(l + Si,o)(l + Sm,o) 
sin(Wp) cos(mt') dt' , 



(41) 



where 9p is the angular position of particle p and N is the 
total number of particles. The integral is over a time period 
of 27r (one orbital period). 8 is the Kronecker delta. The 
other component modes have similar definitions. Four disk 
modes were followed in the simulations: 

(1,0), a measure of the disk eccentricity. 
(3,3), the disk's tidal response to the m = 3 component 
of the binary potential. 
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Figure 7. Grayscale maps of the surface density in simulation 1 at times t = 50.0, 350.0, 400.0, 450.0, 500.0, 550.0 fij^ ^ . The Roche lobe 
and 3 : 1 eccentric inner Lindblad resonance are marked in each frame. A logarithmic scale (to right) is used. 
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Figure 8. Grayscale maps of the surface density in simulation 2 at times t = 50.0, 100.0, 150.0, 200.0, 250.0, 300.0 fij^ ^ . The Roche lobe 
and 3 : 1 eccentric inner Lindblad resonance are marked in each frame. A logarithmic scale (to right) is used. 



12 J. R. Murray 



Figure 9. Fourier mode analysis of disk 1 from time t = 300.00 
to t = 600.00 Qj^ . This covers the interval of strongest growth in 
the disk eccentricity. 



Figure 10. Fourier mode analysis of disk 2 from time t = to 



t 



200.00 fi: 



This covers the interval of strongest growth in 



the disk eccentricity 

(2,3), travelling (two armed spiral) wave launched at the 
Lindblad resonance. This is the third mode in the TRI mech- 
anism. 

(2,2), the disk's tidal response to the m = 2 component 
of the binary potential. Lubow (1992) found that this mode 
cleared material out of the Lindblad resonance and inter- 
fered with the TRI mechanism. 

Figs 9 and 10 plot the various mode strengths as func- 
tions of time for disks 1 and 2 respectively. In both cases the 
graphs cover the period from the disks' initial encounter with 
the resonance up till an eccentric equilibrium is reached. We 



draw attention to two key findings of the analysis. Firstly, 
in both simulations the initial growth in the eccentricity can 
be fit with exponential curves. Secondly the strength of the 
disk's (2, 2) tidal response diminishes over the same time 
interval in which the (1,0) mode gains strength. 

We use Eqn 3 to obtain a theoretical estimate for the 
eccentricity growth rate. The mean radius of the m = 3 
eccentric inner Lindblad resonance is rres = 0.455 d. For 
disk 2, the tidal instability becomes important after about 
t = 50.00 At this time the mean surface density at the 

resonance radius divided by the disk mass 

S(»'re 



M 



1.15 r 



(42) 



The eccentricity is assumed to be independent of radius. The 
theoretical estimate of the eccentricity growth rate is then 
At = 0.097 However, given our assumptions, we feel it is 
probably wise to assume an uncertainty of about 10% and 
take 



At = 0.10 ± 0.01 



(43) 



By simply plotting the (1,0) mode strength on log-linear 
axes and fitting a straight line to the most rapidly increasing 
section of the curve we obtain 



Ai = 0.046 ± 0.002 
for simulation 1, and 
A2 = 0.058 ± 0.002 



(44) 



(45) 



in the second case. The simulation A are unambiguously 
smaller than the theoretical value. However they are only 
reduced from At by a factor 2, which is surprising consider- 
ing that the theory only took account of the m = 3 compo- 
nent of the binary potential and our simulations incorporate 
the full tidal field. In contrast, SPH simulations of Lubow 
(1991b) completed with the full binary potential produced 
only a quadratic growth in the eccentricity. We also see a 
significant difference in eccentricity growth rate between our 
two simulations. In simulation 1, the mass transfer stream 
is feeding mass with low specific angular momentum almost 
directly into the resonance region. It is to be expected there- 
fore that this inhibits the growth of the eccentricity. It does 
not however disrupt the tidal instability mechanism com- 
pletely. 

Lubow (1992) completed a series of simulations with 
only the m = 3, and various fractions of the m = 2 compo- 
nent of the tidal potential included. He found that the tidal 
response generated by the m = 2 component was clearing 
material out of the 3 : 1 resonance and reducing the effec- 
tiveness of the eccentric instability. In our simulations how- 
ever, the (2, 2) mode itself is severely weakened by the disk's 
encounter with the 3 : 1 resonance, and the tidal response's 
capacity to interfere with the eccentricity growth mechanism 
is limited. 

The disks in our simulations were significantly differ- 
ent to those of Lubow, even before the 3 : 1 resonance was 
encountered. We allowed the disk to build up via steady 
mass transfer, and by the time the disk had become reso- 
nant, the inner disk had a density versus radius profile that 
was close to the steady state axisymmetric solution. Lubow 
on the other hand initially laid down an annulus of parti- 
cles on nonintersecting orbits and did not allow any form 
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of mass addition. With viscous and pressure forces present, 
the annulus would spread. Lubow's disk could not approach 
a steady state. For the calculations he performed using only 
the m = 3 component of the tidal potential, the instability 
occurred immediately and grew on a time scale much shorter 
than the spreading time of the ring. However with the full 
potential included, Lubow found he had to wait much longer 
for the eccentric instability to be excited. To overcome that 
problem he gave the annulus an initial, arbitrary eccentric- 
ity. 

The annuli in Lubow's calculations were hotter and less 
viscous than ours, but not dramatically so. He used a sound 
speed c = 0.05 Qb, two and a half times greater than ours. 
His viscous dissipation was characteristic of a disk with a ~ 
0.1. Thus the kinematic viscosity p at resonance would be 
about a factor three less in Lubow's simulations than in ours. 
However Lubow ran an m = 3 simulation with a kinematic 
viscosity increased to be somewhat closer to ours, and found 
no significant difference in his results. 



Table 3. Disk rotation periods obtained in eccentric disk simu- 
lations 



Author 


9 


T 


Comment 


Whitehurst 1988a 


0.15 


1.035 


ly r 2 


Lubow 1991b 


0.10 


1.026 


a ~ 0.1 


Hirose and Osaki 1990 


0.05 


1.044 ± 0.011 


/ = 0.02 




0.1 


1.037 ± 0.021 






0.15 


1.060 ± 0.011 






0.2 


1.069 ± 0.021 






0.25 


1.104 ± 0.028 






0.3 


1.137± 0.024 






0.5 


1.226 ± 0.021 






0.1 


1.043 ± 0.010 


/ = 0.01 




0.125 


1.051 ± 0.010 






0.15 


1.060 ± 0.011 






0.175 


1.081 ± 0.043 






0.2 


1.083 ± 0.025 




Heemskerk 1994 


0.2 


0.97 




Murray 1995 


3 
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1.08 ± 0.02 


Simulation 1 




1.07± 0.02 


Simulation 2 



4.4 Disk precession 

Both disk 1 and disk 2 eventually evolved to a state of ec- 
centric equilibrium. In this state the disks rotated in a retro- 
grade fashion about the white dwarf primary. For simulation 
1 we measured the motion of the disk centre of mass over 
65 cycles to obtain a mean disk period. As a fraction of the 
binary period, ± 3 standard deviations 

^ = 1.08 ± 0.02. (46) 

In simulation 2, over 118 disk cycles 

= 1.07 ± 0.02. (47) 

In neither simulation was there evidence of any secular or 
long term trend in the disk period. The details of the disk 
simulations completed by various authors have already been 
discussed. The precession periods they obtained are listed in 
Table 3. It is apparent that, in addition to any dependance 
of Td upon q, the pressure and viscosity have a significant 
impact upon the precession rate. This comes as a qualitative 
confirmation of Lubow's (1992) identification of the factors 
contributing to the disk's precession. In general, the super- 
hump period is observed to decrease over the course of a 
superoutburst. This is perhaps indicative that the disk prop- 
erties (pressure, temperature, viscosity) in the region of the 
resonance, change significantly as the outburst progresses. 

4.5 Disk dissipation 

Viscous dissipation releases heat in the disk. In an axisym- 
metric thin disk, changes in the radial disk structure occur 
only on the viscous diffusion time scale which is much longer 
than the thermal time scale. In such a disk it can be as- 
sumed that heat is radiated away at the radius at which it 
was generated. In our nonaxisymmetric simulations, gas in 
the outer regions experiences significant changes in condi- 
tions on dynamical time scales which are of similar order to 
thermal time scales. In nonaxisymmetric regions, we cannot 
assume that the disk shines brightest where the viscous dis- 
sipation is greatest. Thus in this section we do not claim to 



Figure 11. Rate of dissipation as a function of time for, the 
entire disk (upper curve), and the outer disk (lower curve) for 
simulation 2. Scaled units are used. The disk is in a precessing 
equilibrium state with an approximately constant mass (see Fig. 
6). 

produce disk luminosity maps, we simply locate the regions 
of significant energy release in these precessing disks. 

In Fig. 11 we plot the rate of viscous dissipation (the 
units are energy time~^) for various regions of disk 2 as a 
function of time. The disk's total dissipation rate (upper- 
most curve) is dominated by a very noisy component from 
the inner disk. In the simulations there are relatively few 
particles near the inner edge where the shear forces (and 
dissipation) are largest. Hence the accretion of even a single 
particle onto the white dwarf makes a large difference to the 
total dissipation rate. 

A second curve, showing the total dissipation rate at 
radii r > 0.05 d is much less noisy. Superposed on a steady 
signal, we see a set of evenly spaced 'humps'. The spacing of 
the humps matches the period of motion of the disk centre of 
mass. The humps are broad, each lasting approximately two 
and a half scaled time units, and have an irregular spiky 
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Figure 12. Grey scale maps of dissipation rate in disk 2 at times t = 1001.00, 1002.00, 1003.00, 1004.00, 1005.00, 1006.00 fi^ ^ . A logarithmic 
scale (shown to right) is used. The inner disk (r < 0.05 d) is not shown. The Roche lobe is marked in each frame. 
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Figure 13. Dissipation rate as a function of time (in days) for 
simulation 2 over the period during whcih the eccentricity growth 
is fastest. The disk encounters the 3 : 1 resonance at approxi- 
mately t = 0.6 days. 



appearance. Recall that in simulation 2 there is no mass 
transfer stream and hence no bright spot. In fact the stream 
impact with the periodically advancing and retreating disk 
edge (seen in simulation 1) produces another series of smaller 
humps that are approximately 180 degrees out of phase with 
those shown in Fig. 11. The first set of humps is simply 
caused by the periodic compression of the eccentric disk as 
it rotates in the white dwarf's Roche lobe. The source is an 
extended region in the outer disk. The best we can say is 
that this is not inconsistent with observational studies. For 
example, Naylor et al. (1987) estimated that the superhump 
light from OY Car came from a region that was optically 
thick with a temperature of approximately 8000K and an 
area comparable to the entire disk. 

Fig. 12 is a sequence of six 'dissipation maps' of disk 2. 
The six frames are separated 1.00 Q,^^ in time and so approx- 
imately cover one full rotation of the disk about the primary. 
The inner portion of the disk r < 0.05 is not shown as the 
large dissipation from this region overwhelms the structure 
in the outer disk. A logarithmic scale is used. 

Frame one shows a relatively compact disk with the 
dissipation dominated by the axisymmetric inner regions. 
Although the outer edge has a lot of structure, there is no 
apparent 'tail' at this stage. 

In frame two (1.00 Qj^^ later) the disk has a large tail 
extending out to three o'clock. This indicates that signifi- 
cant angular momentum and energy has been transferred to 
gas in the outer disk, allowing it to move to larger, more 
eccentric orbits. By the third frame the large tail has moved 
to five o'clock (remember that the motion of individual par- 
ticles is anticlockwise). In frame four the tail has moved 
further clockwise to about seven o'clock. We note the ap- 
pearance of a narrow dense, high dissipation region along 
the trailing edge of the tail. In the arc from about seven 
to nine o'clock, the outer disk experiences a negative tidal 
torque (see e.g. Fig. 3.6 Hirose 1991). In other words, in this 
region angular momentum is extracted from the outer disk 
and returned to the binary orbit. Thus, as the gas is braked 



by this tidal torque, we see a compression at the trailing 
edge of the tail. 

By frame 5 the tail has moved closer to eight o'clock and 
the compression at the trailing edge has intensified. Consid- 
erable amounts of kinetic energy are being released as heat 
here. In between frames four and five the disk's dissipation 
rate has increased. Frame 6 corresponds to maximum dissi- 
pation rate. We see the tail generated over the last precession 
cycle in the process of being destroyed. In particular we see 
one of the dense knots that form in the trailing edge of the 
tail, falling into the axisymmetric inner disk, releasing large 
amounts of heat. In fact the knot shown here is responsi- 
ble for the second largest dissipation spike in this particular 
hump. 

Over the course of one rotation, the disk develops a 
considerable eccentric tail. This tail is then destroyed by 
negative tidal torques, releasing angular momentum to the 
binary orbit and energy as heat. 

O'Donoghue (1990) attempted to use the eclipse map- 
ping technique of Home (1985) to obtain two dimensional 
maps of the superhump source. The original technique was 
only used to map sources that varied on time scales much 
longer than the eclipse time. In order to apply the tech- 
nique to sources that are variable on shorter time scales, 
additional constraints upon the problem are needed. Hence 
O'Donoghue assumed that the superhump source had a fixed 
position in the binary frame, and that the entire superhump 
source varied in brightness coherently. With these assump- 
tions, O'Donoghue considered four eclipsed superhumps in 
detail. In each case he found the source to be located on 
the edge of the disk in the approximate direction of the 
secondary star (i.e. towards 9 o'clock if we were looking at 
Fig. 12). In three eclipses the source could be resolved into 
three components; a compact source at the disk edge lying 
almost on the disk axis, and two flanking extended sources 
that followed the disk edge in either direction. The fourth 
eclipse showed a single source located near where the mass 
transfer stream would be expected to strike the disk edge. 

There are significant differences between our two di- 
mensional simulations, and the superhump eclipse maps of 
O'Donoghue. Whereas he has determined that the super- 
hump source is located close to the binary axis, we see the 
dissipation occurring downstream (anticlockwise) of this lo- 
cation. In the simulations, the sources of excess dissipation 
are not stationary but move with the gas flow. Neither are 
they a single coherent source. It would be interesting to see 
how those two assumptions affected O'Donoghue's results. 
On the other hand, we have not treated the thermodynam- 
ics and radiation physics at all in these simulations, and 
we do not portray our dissipation maps (Fig. 12) as being 
equivalent to luminosity maps. 

As a coda to this section, we show that the tidal reso- 
nance instability can generate a signal on a time scale con- 
sistent with observations. Superhumps first appear in the 
lightcurve of an SU UMa system within a day of it going 
into superoutburst. In Fig. 13 we show the dissipation rate 
in the outer regions (r > 0.05) of disk 2 as a function of time 
in days (recall that we chose system parameters appropriate 
to OY Car in superoutburst). The disk first encountered the 
resonance at t ~ 0.5 days. The graph clearly shows a peri- 
odic signal growing to full strength within the next day or 
so. As we have said before, these are very simple simulations 
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that concentrate on the dynamics of the system. The details 
of the disk thermodynamics and radiation have been largely 
ignored. That a significant periodic signal does appear on 
a similar time scale to superhump observations is not proof 
that the precessing disk model is correct. If however, the 
periodicity in the dissipation rate had only emerged after a 
much longer time, then we would suspected that the tidal in- 
stability mechanism was too weak to generate superhumps. 



5 SUMMARY 

Particle simulations of accretion disks in close binary sys- 
tems have been reviewed. Particle methods are well suited 
to modelling highly viscous disks in gravitational potentials 
with little symmetry. A Smoothed Particle Hydrodynamics 
(SPH) scheme specifically for thin disk problems was devel- 
oped. We described the use of an artificial viscosity term 
to provide viscous dissipation and then derived a theoreti- 
cal estimate for an equivalent shear viscosity. Axisymmetric 
and nonaxisymmetric tests of the code confirmed the theo- 
retical result to within 10%. This is an improvement upon 
earlier particle schemes for which the dissipation could only 
be determined empirically. 

Two disk simulations for a mass ratio q = j^, differ- 
ing only in the method of mass addition to the disk, were 
completed. We find that the addition of material with a low 
specific angular momentum to the outer disk, via the mass 
transfer stream from Li , circularises disk orbits and inhibits 
the disk's radial expansion. 

Both our disks developed an eccentricity in a fashion 
consistent with the tidal resonance mechanism detailed in 
Lubow (1991a). The initial growth of the eccentric mode 
strength was exponential with a rate approximately half the 
predicted value. Mass transfer via the Li stream reduced 
the growth rate slightly. Lubow (1992) suspected that the 
disk's tidal response to the m = 2 component of the binary 
potential reduced the strength of the eccentric instability by 
removing material from the resonance region. In our simu- 
lations, the m = 2 tidal response weakened considerably as 
the disk became more eccentric. 

In the binary frame, the eccentric disks rotated in a 
retrograde fashion with a period slightly exceeding the bi- 
nary's. The rotation gave rise to a periodic excess in the 
rate at which energy is viscously dissipated in the disks. 
In the precessing disk model, this excess is ultimately the 
source of the superhump luminosity. We find that as the 
resonant disk rotates, a large eccentric tail develops. When 
the tail moves through a region where the tidal torques are 
negative it collapses in toward the centre of the disk, releas- 
ing large amounts of heat. The source of the dissipation is 
found to share general properties with superhump observa- 
tions however the details cannot be immediately reconciled 
with the eclipse maps of O'Donoghue (1990). Our results are 
certainly inconsistent with the assumptions used to gener- 
ate those maps. It is hoped that these simulations will be 
useful in determining new constraints for the deconvolution 
of superhump eclipses. 

Finally we find that the periodic component of the dis- 
sipation curve developed on a time scale that was consistent 
with the first appearance of superhumps in a superoutburst. 
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